# -*- coding: utf-8 -*-
"""
@authors: coude

This library is designed to create the figures for the FIELDMAPS Data Release paper. 
    
"""

# REMINDER: Python arrays are inverted [y,x] relative to IDL [x,y]

# =============================================================================
# Package dependencies
# =============================================================================

from aplpy import FITSFigure
from astropy.io import fits
from astropy.wcs import WCS
from lmfit import Model
import matplotlib.patches as patches
import matplotlib.pyplot as plt
import math
import numpy as np
from scipy import stats, optimize


# =============================================================================
# Flux estimate at a given frequency
# =============================================================================

def flux_calculation(nu,T,N):
    # Variables
    # nu : frequency in Hz
    # T : temperature in K
    # Column density in cm**-2
    # Constants
    h = 6.626*10**(-34.0) # Planck constant in J s
    c = 2.998*10**8.0  # Speed of light in m/s
    k = 1.381*10**(-23.0) # Boltzmann constant in J/K
    mu = 2.8 # Mean molecular weigth
    beta = 1.75 # Spectral index of emissivity
    mh = 1.6738*10**(-27.0)  # Mass hydrogen in kg
    Rgd = 100 # Gas to dust ratio

    # Converting density to m**-2
    NH = N*10000.0

    A1=2*h/c**2
    A2 = h/k
    exp1 = A2*nu/T
    exp2 = np.exp(exp1)
    Bv = A1*nu**3*(exp2-1)**-1 # Planck law in W sr^-1 m^-2 Hz^-1

    Kv = 4.0 * (nu / (505*10**9))**beta * 0.1 # Opacity law m^2 kg^-1 (Wang 2014)
    Tau = mu * mh * Kv * NH / Rgd # Opacity (Wang 2014), mu from Kauffmann 2008

    Iv = Bv * (1-np.exp(-Tau)) # Flux in W sr^-1 m^-2 Hz^-1
    IvJy = 1000.0 * 10**26 * Iv * (4.255*10**10)**-1 # Flux in mJy arcsec^-2

    # Printing the result
    #print('The requested flux is: '+str(IvJy) + ' Jy arcsec**-2')
    return IvJy

# =============================================================================
# Polarization map on temperature
# =============================================================================

def CreateFluxMap214(Nref,Tref):
    nu = 1.4*10**12 # 214 µm in Hz
    # Measuring the size of the new array
    shapeX = Nref.data.shape[0]
    shapeY = Nref.data.shape[1]
    # Creating new array
    fake_flux = np.zeros((shapeX,shapeY))
    # Calculating estimated flux for each pixel
    for i in range(0,shapeX):
        for j in range(0,shapeY):
            fake_flux[i,j] = flux_calculation(nu,Tref.data[i,j],Nref.data[i,j])
    # Create FITS file
    FluxMap = fits.PrimaryHDU(data=fake_flux, header=Nref.header)
    return FluxMap

# =============================================================================
# Magnetic field map on column density
# =============================================================================

def BvNmap(Nref, Iref, Pref, ScaleLength, Region, FigSize, Ilevels, Nscale, IdI=10.0, PdP=3.0, Pmax=30.0,
           StepVec=3, ScaleVec=2.0, Colorbar='top', BeamPos='bottom left', BeamPad=1.0, Scalebar='top right', Linewidth=2.0, HAWCBeam=0.00506, 
           FilLabel=None, FilLabelX=0.5, FilLabelY=0.95):
    # Creating a copy of the column density FITS container to have a cleaner color scale label
    Narray = Nref.data.copy()
    Narray = Narray*10**-22
    Nfits = fits.PrimaryHDU(data=Narray, header=Nref.header)
    # Plotting the COLUMN DENSITY MAGNETIC FIELD MAP
    BvN_plot = FITSFigure(Nfits,figsize=(FigSize[0], FigSize[1]))
    BvN_plot.tick_labels.set_xformat('dd.d')
    BvN_plot.tick_labels.set_yformat('d.dd')
    BvN_plot.show_colorscale(cmap='Greys', vmin=Nscale[0]*10**-22, vmax=Nscale[1]*10**-22)
    BvN_plot.add_colorbar(pad=0.25)
    BvN_plot.colorbar.set_location(Colorbar)
    BvN_plot.colorbar.set_axis_label_text(r'$\it{Herschel}$ Column Density $N_{H_2}$ ($10^{22}$ cm$^{-2}$)')
    BvN_plot.add_scalebar(ScaleLength, '1 pc', corner=Scalebar, frame=True) # Scalebar equivalent for 1 pc 
    # Adding beams
    # Herschel
    BvN_plot.add_beam(facecolor='green', edgecolor='green', zorder=1,
                        linewidth=2, pad=1.0, corner=BeamPos) # Herschel Beam
    BvN_plot.beam.set_major(0.01011) # Update major axis for Herschel beam (36.4'', see Aniano et al. 2011)
    BvN_plot.beam.set_minor(0.01011) # Update minor axis for Herschel beam 
        # HAWC+
    BvN_plot.add_beam(facecolor='red', edgecolor='red',
                    linewidth=2, pad=BeamPad, corner=BeamPos) # HAWC+ Beam
    BvN_plot.beam[1].set_major(HAWCBeam) # Update major axis for HAWC+ Beam
    BvN_plot.beam[1].set_minor(HAWCBeam) # Update minor axis for HAWC+ Beam

    # Adding total intensity contour
    BvN_plot.show_contour(data=Iref, colors='white', levels=Ilevels)

    # Centering
    BvN_plot.recenter(x=Region[0], y=Region[1], width=Region[2], height=Region[3])

    # Creating a new fits object for APLpy's FITSFigure method to recognize
    # for the polarization fraction P
    pref = fits.PrimaryHDU(data=Pref[8].data, header=Pref[0].header)
    pmap = pref.copy()
    # for the magnetic field angle B
    oref = fits.PrimaryHDU(data=Pref[11].data, header=Pref[0].header)
    omap = oref.copy()

    # Creating a mask to hide polarization vectors with low signal-to-noise ratios
    imask_01 = np.where(Pref[0].data/Pref[1].data < IdI) # Total intensity SNR threshold
    imask_02 = np.where(Pref[0].data < 0) # Total intensity positive threshold
    pmask = np.where(Pref[8].data/Pref[9].data < PdP) # Polarization SNR threshold
    pmaxmask = np.where(Pref[8].data > Pmax) # Polarization SNR threshold
    # Forcing the polarization vectors to share the same amplitude
    pmap.data[np.where(Pref[8].data > 0.0)] = 1.0

    # Masking all the indices for  which the selection criteria failed
    pmap.data[imask_01] = np.nan
    pmap.data[imask_02] = np.nan
    pmap.data[pmask] = np.nan
    pmap.data[pmaxmask] = np.nan

    # Plotting the polarization vectors
    BvN_plot.show_vectors(pmap, omap, scale=ScaleVec, step=StepVec, color='red', zorder=3, alpha=0.8, linewidth=Linewidth)

    # Adding label for the filament
    if FilLabel != None:
        BvN_plot.add_label(FilLabelX, FilLabelY, FilLabel, relative=True, size=13)


    return BvN_plot

# =============================================================================
# Adding Planck vectors to magnetic field map
# =============================================================================

def BvN_Add_Planck(BvN_plot, Planck, StepVec=3, ScaleVec=2.0, Color='blue', Linewidth=2.0):
    # Creating fits containers
    pmap = Planck[4].copy()
    pmap.data[:,:] = 1.0 # Force length of vectors to be identical
    bmap = Planck[6].copy()
    # Plotting the polarization vectors
    BvN_plot.show_vectors(pmap, bmap, scale=ScaleVec, step=StepVec, color=Color, linewidth=Linewidth)
    # Returning the updated plot (not required, only for clarity)
    return BvN_plot

# =============================================================================
# Polarization map on Stokes I total intensity
# =============================================================================

def PvImap(Iref, Pref, Region, FigSize, Ilevels, Iscale, IdI=10.0, PdP=3.0, Pmax=30.0,
           StepVec=3, ScaleVec=0.8, Colorbar='top', BeamPos='bottom left', Scalebar='top right', HAWCBeam=0.00506,
           FilLabel=None, FilLabelX=0.5, FilLabelY=0.95):
    # Plotting the INTENSITY POLARIZATION MAP
    PvI_plot = FITSFigure(Iref,figsize=(FigSize[0], FigSize[1]))
    PvI_plot.tick_labels.set_xformat('dd.d')
    PvI_plot.tick_labels.set_yformat('d.dd')
    PvI_plot.show_colorscale(cmap='Greys', vmin=Iscale[0], vmax=Iscale[1])
    PvI_plot.add_colorbar(pad=0.25)
    PvI_plot.colorbar.set_location(Colorbar)
    PvI_plot.colorbar.set_axis_label_text('HAWC+ Stokes $I$ Intensity (mJy arcsec$^{-2}$)')
    PvI_plot.add_beam(facecolor='red', edgecolor='red',
                        linewidth=2, pad=1, corner=BeamPos) # HAWC+ Beam
    PvI_plot.beam.set_major(HAWCBeam) # Update major axis for HAWC+ Beam
    PvI_plot.beam.set_minor(HAWCBeam) # Update minor axis for HAWC+ Beam
    # Adding the vector scale bar
    vectscale = ScaleVec * Iref.header['CDELT2'] # AplPy calculates vector scale with pixel units
    PvI_plot.add_scalebar(5 * vectscale, "P = 5%",corner=Scalebar,frame=True) # Scalebar equivalent for 1 pc 

    # Adding total intensity contour
    PvI_plot.show_contour(data=Iref, colors='white', levels=Ilevels)
    # Centering
    PvI_plot.recenter(x=Region[0], y=Region[1], width=Region[2], height=Region[3])

    # Creating a new fits object for APLpy's FITSFigure method to recognize
    # for the polarization fraction P
    pref = fits.PrimaryHDU(data=Pref[8].data, header=Pref[0].header)
    pmap = pref.copy()
    # for the magnetic field angle B
    oref = fits.PrimaryHDU(data=Pref[10].data, header=Pref[0].header)
    omap = oref.copy()
    # Creating a mask to hide polarization vectors with low signal-to-noise ratios
    imask_01 = np.where(Pref[0].data/Pref[1].data < IdI) # Total intensity SNR threshold
    imask_02 = np.where(Pref[0].data < 0) # Total intensity positive threshold
    pmask = np.where(Pref[8].data/Pref[9].data < PdP) # Polarization SNR threshold
    pmaxmask = np.where(Pref[8].data > Pmax) # Polarization SNR threshold
    # Masking all the indices for  which the selection criteria failed
    pmap.data[imask_01] = np.nan
    pmap.data[imask_02] = np.nan
    pmap.data[pmask] = np.nan
    pmap.data[pmaxmask] = np.nan

    # Plotting the polarization vectors
    PvI_plot.show_vectors(pmap, omap, scale=ScaleVec, step=StepVec, color='red')

    # Adding label for the filament
    if FilLabel != None:
        PvI_plot.add_label(FilLabelX, FilLabelY, FilLabel, relative=True, size=13)

    return PvI_plot

# =============================================================================
# Polarization map on column density
# =============================================================================

def PvNmap(Nref, Iref, Pref, Region, FigSize, Ilevels, Nscale, IdI=10.0, PdP=3.0, Pmax=30.0,
           StepVec=3, ScaleVec=0.8, Colorbar='top', BeamPos='bottom left', HAWCBeam=0.00506):
    # Plotting the COLUMN DENSITY POLARIZATION MAP
    PvN_plot = FITSFigure(Nref,figsize=(FigSize[0], FigSize[1]))
    PvN_plot.tick_labels.set_xformat('dd.d')
    PvN_plot.tick_labels.set_yformat('d.dd')
    PvN_plot.show_colorscale(cmap='Greys', vmin=Nscale[0], vmax=Nscale[1])
    PvN_plot.add_colorbar(pad=0.25)
    PvN_plot.colorbar.set_location(Colorbar)
    PvN_plot.colorbar.set_axis_label_text('Column Density $N_{H_2}$ (cm$^{-2}$)')
    # Adding the vector scale bar
    vectscale = ScaleVec * Iref.header['CDELT2'] # AplPy calculates vector scale with pixel units
    PvN_plot.add_scalebar(5 * vectscale, "P = 5%",corner='top right',frame=True) # Scalebar equivalent for 1 pc 

    # Adding beams
    # Herschel
    PvN_plot.add_beam(facecolor='green', edgecolor='green', zorder=1,
                        linewidth=2, pad=1, corner=BeamPos) # Herschel Beam
    PvN_plot.beam.set_major(0.01011) # Update major axis for Herschel beam (36.4'', see Aniano et al. 2011)
    PvN_plot.beam.set_minor(0.01011) # Update minor axis for Herschel beam 
    # HAWC+
    PvN_plot.add_beam(facecolor='red', edgecolor='red',
                    linewidth=2, pad=1, corner=BeamPos) # HAWC+ Beam
    PvN_plot.beam[1].set_major(HAWCBeam) # Update major axis for HAWC+ Beam
    PvN_plot.beam[1].set_minor(HAWCBeam) # Update minor axis for HAWC+ Beam

    # Adding total intensity contour
    PvN_plot.show_contour(data=Iref, colors='white', levels=Ilevels)
    # Centering
    PvN_plot.recenter(x=Region[0], y=Region[1], width=Region[2], height=Region[3])

    # Creating a new fits object for APLpy's FITSFigure method to recognize
    # for the polarization fraction P
    pref = fits.PrimaryHDU(data=Pref[8].data, header=Pref[0].header)
    pmap = pref.copy()
    # for the magnetic field angle B
    oref = fits.PrimaryHDU(data=Pref[10].data, header=Pref[0].header)
    omap = oref.copy()
    # Creating a mask to hide polarization vectors with low signal-to-noise ratios
    imask_01 = np.where(Pref[0].data/Pref[1].data < IdI) # Total intensity SNR threshold
    imask_02 = np.where(Pref[0].data < 0) # Total intensity positive threshold
    pmask = np.where(Pref[8].data/Pref[9].data < PdP) # Polarization SNR threshold
    pmaxmask = np.where(Pref[8].data > Pmax) # Polarization SNR threshold
    # Masking all the indices for  which the selection criteria failed
    pmap.data[imask_01] = np.nan
    pmap.data[imask_02] = np.nan
    pmap.data[pmask] = np.nan
    pmap.data[pmaxmask] = np.nan

    # Plotting the polarization vectors
    PvN_plot.show_vectors(pmap, omap, scale=ScaleVec, step=StepVec, color='red')

    return PvN_plot

# =============================================================================
# Polarization map on temperature
# =============================================================================

def PvTmap(Tref, Iref, Pref, Region, FigSize, Ilevels, Tscale, IdI=10.0, PdP=3.0, Pmax=30.0,
           StepVec=3, ScaleVec=0.8, Colorbar='top', BeamPos='bottom left', HAWCBeam=0.00506, Scalebar='top right'):
    # Plotting the TEMPERATURE POLARIZATION MAP
    PvT_plot = FITSFigure(Tref,figsize=(FigSize[0], FigSize[1]))
    PvT_plot.tick_labels.set_xformat('dd.d')
    PvT_plot.tick_labels.set_yformat('d.dd')
    PvT_plot.show_colorscale(cmap='rainbow', vmin=Tscale[0], vmax=Tscale[1])
    PvT_plot.add_colorbar(pad=0.25)
    PvT_plot.colorbar.set_location(Colorbar)
    PvT_plot.colorbar.set_axis_label_text('Dust Temperature $T_d$ (K)')
    # Adding the vector scale bar
    vectscale = ScaleVec * Iref.header['CDELT2'] # AplPy calculates vector scale with pixel units
    PvT_plot.add_scalebar(5 * vectscale, "P = 5%",corner=Scalebar,frame=True) # Scalebar equivalent for 1 pc 

    # Adding total intensity contour
    PvT_plot.show_contour(data=Iref, colors='white', levels=Ilevels)
    # Centering
    PvT_plot.recenter(x=Region[0], y=Region[1], width=Region[2], height=Region[3])

    # Adding beams
    # Herschel
    PvT_plot.add_beam(facecolor='red', edgecolor='red', zorder=1,
                        linewidth=2, pad=1, corner=BeamPos) # Herschel Beam
    PvT_plot.beam.set_major(0.01011) # Update major axis for Herschel beam (36.4'', see Aniano et al. 2011)
    PvT_plot.beam.set_minor(0.01011) # Update minor axis for Herschel beam 
        # HAWC+
    PvT_plot.add_beam(facecolor='black', edgecolor='black',
                    linewidth=2, pad=1, corner=BeamPos) # HAWC+ Beam
    PvT_plot.beam[1].set_major(HAWCBeam) # Update major axis for HAWC+ Beam
    PvT_plot.beam[1].set_minor(HAWCBeam) # Update minor axis for HAWC+ Beam

    # Creating a new fits object for APLpy's FITSFigure method to recognize
    # for the polarization fraction P
    pref = fits.PrimaryHDU(data=Pref[8].data, header=Pref[0].header)
    pmap = pref.copy()
    # for the magnetic field angle B
    oref = fits.PrimaryHDU(data=Pref[10].data, header=Pref[0].header)
    omap = oref.copy()
    # Creating a mask to hide polarization vectors with low signal-to-noise ratios
    imask_01 = np.where(Pref[0].data/Pref[1].data < IdI) # Total intensity SNR threshold
    imask_02 = np.where(Pref[0].data < 0) # Total intensity positive threshold
    pmask = np.where(Pref[8].data/Pref[9].data < PdP) # Polarization SNR threshold
    pmaxmask = np.where(Pref[8].data > Pmax) # Polarization SNR threshold
    # Masking all the indices for  which the selection criteria failed
    pmap.data[imask_01] = np.nan
    pmap.data[imask_02] = np.nan
    pmap.data[pmask] = np.nan
    pmap.data[pmaxmask] = np.nan

    # Plotting the polarization vectors
    PvT_plot.show_vectors(pmap, omap, scale=ScaleVec, step=StepVec, color='black')

    return PvT_plot

# =============================================================================
# Creating the polarization vector catalog
# =============================================================================

def MakeCats(Nref, Tref, Iref, Pref, Iest, CatMask=None, IdI=10.0, PdP=3.0, Pmax=30.0):
    Ndata = Nref.copy()
    Tdata = Tref.copy()
    Idata = Iref.copy()
    Ifake = Iest.copy()
    Pdata = Pref.copy()

    # Creating lists
    CatI = [] # Stokes I
    CatIest = [] # Stokes I as estimated from Herschel
    CatdI =[] # Uncertainty dI for Stokes I
    CatQ =[] # Stokes Q
    CatdQ =[] # Uncertainty dQ for Stokes Q
    CatU =[] # Stokes U
    CatdU =[] # Uncertainty dU for Stokes U
    CatP = [] # Polarization fraction P (debiased)
    CatdP =[] # Uncertainty dP for polarization fraction P
    CatO = [] # Polarization angle O
    CatB = [] # Rotated polarization angle B (+90 degrees)
    CatdO =[] # Uncertainty dO for polarization angles
    CatPI = [] # Polarized intensity PI (debiased)
    CatdPI =[] # Uncertainty dPI for polarized intensity PI
    CatN = [] # Column density
    CatT =[] # Dust temperature

    fluxconv = 1000.0*(3600.0*Idata.header['CDELT2'])**-2.0 # Conversion factor from Jy/pixel

    # Creating a new fits object for APLpy's FITSFigure method to recognize
    # for the polarization fraction P
    pref = fits.PrimaryHDU(data=Pdata[8].data, header=Pdata[0].header)
    pmap = pref.copy()
    # for the magnetic field angle B
    oref = fits.PrimaryHDU(data=Pdata[10].data, header=Pdata[0].header)
    omap = oref.copy()
    # Creating a mask to hide polarization vectors with low signal-to-noise ratios
    imask_01 = np.where(Pdata[0].data/Pdata[1].data < IdI) # Total intensity SNR threshold
    imask_02 = np.where(Pdata[0].data < 0) # Total intensity positive threshold
    pmask = np.where(Pdata[8].data/Pdata[9].data < PdP) # Polarization SNR threshold
    pmaxmask = np.where(Pdata[8].data > Pmax) # Polarization SNR threshold
    # Masking all the indices for  which the selection criteria failed
    pmap.data[imask_01] = np.nan
    pmap.data[imask_02] = np.nan
    pmap.data[pmask] = np.nan
    pmap.data[pmaxmask] = np.nan

    # If region mask is provided
    if CatMask != None:
        Maskdata = CatMask.copy() # Mask for the catalog
        cmask = np.where(Maskdata.data != 1.0) # Catalog mask
        pmap.data[cmask] = np.nan

    # Measuring the size of the new array
    shapeX = Ndata.data.shape[0]
    shapeY = Ndata.data.shape[1]

    for i in range (0,shapeX):
        for j in range (0,shapeY):
            if pmap.data[i,j] > 0.0: # Iterate only on pixels where the selection criteria are respected
                CatI.append(Idata.data[i,j])
                CatIest.append(Ifake.data[i,j])
                CatdI.append(Pdata[1].data[i,j]* fluxconv)  # Conversion factor from Jy/pixel
                CatQ.append(Pdata[2].data[i,j]* fluxconv)
                CatdQ.append(Pdata[3].data[i,j]* fluxconv)
                CatU.append(Pdata[4].data[i,j]* fluxconv)
                CatdU.append(Pdata[5].data[i,j]* fluxconv)
                CatP.append(pmap.data[i,j])
                CatdP.append(Pdata[9].data[i,j])
                CatO.append(Pdata[10].data[i,j])
                CatB.append(Pdata[11].data[i,j])
                CatdO.append(Pdata[12].data[i,j])
                CatPI.append(Pdata[15].data[i,j]* fluxconv)
                CatdPI.append(Pdata[14].data[i,j]* fluxconv)
                CatN.append(Ndata.data[i,j])
                CatT.append(Tdata.data[i,j])

    # Transforming lists into numpy arrays
    CatI = np.array(CatI)
    CatIest = np.array(CatIest)
    CatdI = np.array(CatdI)
    CatQ = np.array(CatQ)
    CatdQ = np.array(CatdQ)
    CatU = np.array(CatU)
    CatdU = np.array(CatdU)
    CatP = np.array(CatP)
    CatdP = np.array(CatdP)
    CatO = np.array(CatO)
    CatB = np.array(CatB)
    CatdO = np.array(CatdO)
    CatPI = np.array(CatPI)
    CatdPI = np.array(CatdPI)
    CatN = np.array(CatN)
    CatT = np.array(CatT)


    # Returning all the catalogs 
    return CatI, CatIest, CatdI, CatQ, CatdQ, CatU, CatdU, CatP, CatdP, CatO, CatB,CatdO, CatPI, CatdPI, CatN, CatT

# =============================================================================
# Creating a mask for a region of the map
# =============================================================================
def MaskRegion(FitsFile, CatRegion, angle=0):
       # Reading region information from CatRegion
       center_l = CatRegion[0]
       center_b = CatRegion[1]
       width = CatRegion[2]
       height = CatRegion[3]
       
       # Creating the empty array for the final mask
       maskref = fits.PrimaryHDU(data=FitsFile[0].data, header=FitsFile[0].header)
       mask = maskref.copy()
       # Size of the map
       MaskSize = np.shape(FitsFile[0].data) 

       # Creating WCS object
       MaskWCS = WCS(FitsFile[0].header)          
       # Identifying pixels to be included in the mask
       for i in range (0,MaskSize[0]): # Latitude
              for j in range (0,MaskSize[1]): # Longitude
                     # Determining the current Galactic position in the array
                     Maskl, Maskb = (MaskWCS.array_index_to_world_values(i, j))
                     # Distance to central pixel
                     delta_l = (Maskl - 
                                   center_l) * math.cos(
                                   Maskb*math.pi/180.0)
                     delta_b = Maskb - center_b

                     # Coordinate conversion
                     x1 = delta_l * -1.0
                     y1 = delta_b
                     ang = angle*math.pi/180.0
                     # Measuring angular distances
                     delta_x2 = abs(x1*math.cos(ang)+y1*math.sin(ang))
                     delta_y2 = abs(-x1*math.sin(ang)+y1*math.cos(ang))
                     # Is the pixel in the mask? If yes, set to 1.0, if no, set to NaN
                     if delta_x2 < 0.5*width and delta_y2 < 0.5*height:
                            mask.data[i,j] = 1.0
                     else:
                            mask.data[i,j] = np.nan
       return mask

# =============================================================================
# Creating a mask to downsample the data for the vector catalogs
# =============================================================================
def DownsampleVectors(FitsFile, Step=1):
    # Creating the empty array for the final mask
    maskref = fits.PrimaryHDU(data=FitsFile[0].data, header=FitsFile[0].header)
    mask = maskref.copy()
    # Set all elements to NaN
    mask.data.fill(np.nan)
    # Size of the map
    MaskSize = np.shape(FitsFile[0].data) 
    for i in range (0,MaskSize[0],Step): # Latitude
            for j in range (0,MaskSize[1],Step): # Longitude
                mask.data[i,j] = 1.0 # Setting every pixel at a given number of steps to 1.0
    return mask

# =============================================================================
# Polarization fraction as a function of total intensity
# =============================================================================
def PvI(CatP, CatdP, CatI, Iscale=None, Pscale=None, showfit='false', weighted='true', errorbars='false', FigSize=[5.12,3.84], Source=None):  
    PvI_plot = plt.figure(figsize=(FigSize[0],FigSize[1])) # Creating the figure object
    plt.xlabel('Stokes Total Intensity $I$ (mJy arcsec$^{-2}$)')
    plt.ylabel("Polarization Fraction $P$ (%)")
    # Plots with and without errorbars 
    if errorbars=='true':
        markers, caps, bars = plt.errorbar(CatI, CatP, yerr=CatdP, fmt ='.', color='black', ecolor='gray', markersize=1, mfc='none',elinewidth=0.9)
        [bar.set_alpha(0.25) for bar in bars]
    else:
        plt.plot(CatI, CatP, '.', color='black', markersize=1, mfc='none')         
    # Switching the total intensity in log scale
    plt.xscale('log')
    axes = plt.gca() # Calling the axes object of the figure
    axes.xaxis.set_ticks_position('both') # Adding ticks to each side
    axes.yaxis.set_ticks_position('both') # Adding ticks to each side
    plt.tight_layout() # Using all available space in the plot window

    # Checking the automated limits for I and P
    xmin, xmax = plt.xlim()
    ymin, ymax = plt.ylim()
    # Setting mininum and maximum values of the I plot
    if Iscale==None:
        Iscale = np.empty(2)
        Iscale[0] = xmin
        Iscale[1] = xmax
    if Pscale==None:
        Pscale = np.empty(2)
        Pscale[0] = ymin
        Pscale[1] = ymax
    # Setting the minimum and maximum values
    plt.xlim(Iscale[0],Iscale[1])
    plt.ylim(Pscale[0],Pscale[1])


    if showfit=='true':
        # Create the array with ln(I), ln(P), and d(ln(P))
        # The numpy package uses log(x) as the natural logarithm ln
        length = CatP.size
        PvIarray = np.empty([length, 3])
        PvIarray[:,0]=np.log(CatI)
        PvIarray[:,1]=np.log(CatP)
        PvIarray[:,2]= CatdP/CatP # Uncertainties in log scale
        # Defining linear function to be fitted by curve_fit
        def lnPvlnI(x, a, b):
            return a*x + b
        # Defining power-law function to be fitted by curve_fit
        def PvItest(x, a, b):
            return b*x**a
        # Fitting power law to data
        if weighted=='true':
            PvIpow, PvIcov = optimize.curve_fit(lnPvlnI, PvIarray[:,0], PvIarray[:,1], sigma=PvIarray[:,2])
            Testpow, Testcov = optimize.curve_fit(PvItest, CatI, CatP, sigma=CatdP)
            # Printing test result
            print('===')
            print('Curve fit test')
            print('Log scale fit: ('+str(PvIpow[0])+','+str(math.exp(PvIpow[1]))+ ')')
            print('Direct power law fit: ('+str(Testpow[0])+','+str(Testpow[1])+ ')')
            print('===')
            print('')
        else:
            PvIpow, PvIcov = optimize.curve_fit(lnPvlnI, PvIarray[:,0], PvIarray[:,1])
            Testpow, Testcov = optimize.curve_fit(PvItest, CatI, CatP)
            # Printing test result
            print('===')
            print('Curve fit test')
            print('Log scale fit: ('+str(PvIpow[0])+','+str(math.exp(PvIpow[1]))+ ')')
            print('Direct power law fit: ('+str(Testpow[0])+','+str(Testpow[1])+ ')')
            print('===')
            print('')
        # Uncertainties on the fit
        PvIerr = np.sqrt(np.diag(PvIcov))
        # Converting the coefficient back to its non-logarithmic value
        PvIcoeff = math.exp(PvIpow[1])
        PvIdcoeff = math.exp(PvIpow[1]) * PvIerr[1]
        # Printing the results of the fit
        print('Power law index: '+str(PvIpow[0])+' ± '+str(PvIerr[0]))
        print('Coefficient: '+str(PvIcoeff)+' ± '+str(PvIdcoeff))
        # Calculating the reduced Chi-Squared value
        if weighted=='false':
            FakeData = lnPvlnI(PvIarray[:,0],PvIpow[0],PvIpow[1])
            Chisq = np.sum((PvIarray[:,0]-FakeData)**2)
            print('Chi-Squared: ' + str(Chisq))
            print('Number of elements: ' + str(PvIarray[:,0].size))
            RedChisq = Chisq/(PvIarray[:,0].size-2) # Reduced Chi-Squared is number of elements minus number of fitted parameters
            print('Reduced Chi-Squared: ' +str(RedChisq))

        if Source != None:
            PvI_plot.text(0.92, 0.85, Source, verticalalignment='bottom', horizontalalignment='right', fontweight='bold')

    # Adding power law fit
    if showfit=='true':
        def PvIfunc(x, A, B):
            return A*x**B
            
        # Creating x values
        PvIfuncX = np.arange(Iscale[0], Iscale[1], 0.1)
        # Calculating the y values
        PvIfuncY = PvIfunc(PvIfuncX,PvIcoeff,PvIpow[0])
        plt.plot(PvIfuncX, PvIfuncY, 'red')
        # Rounding up the fit error
        RoundErrorPvI = round_decimals_up(PvIerr[0],2)
        RoundErrorCoeff = round_decimals_up(PvIdcoeff,1)

        # # Add cases where B=0.5 and B=1.0
        # # Creating x values
        # PvIfuncX_05 = np.arange(Iscale[0], Iscale[1], 0.1)
        # # Calculating the y values
        # PvIfuncY_05 = PvIfunc(PvIfuncX,PvIcoeff,-0.5)
        # # Creating x values
        # PvIfuncX_10 = np.arange(Iscale[0], Iscale[1], 0.1)
        # # Calculating the y values
        # PvIfuncY_10 = PvIfunc(PvIfuncX,PvIcoeff,-1.0)
        # # Plotting the -0.5 index curve
        # plt.plot(PvIfuncX_05, PvIfuncY_05, 'red', linestyle='-.')
        # # Plotting the -1.0 index curve
        # plt.plot(PvIfuncX_10, PvIfuncY_10, 'red', linestyle='--')

        # Adding the fit results as a label
        if weighted=='false':
            Coefficients = '\n'.join((r'$\alpha = {0:.2f} \pm {1:.2f}$'.format(PvIpow[0],RoundErrorPvI),
                                    r'$A = {0:.1f} \pm {1:.1f}$'.format(PvIcoeff,RoundErrorCoeff),
                                    r'$\chi^2_r = {0:.2f}$'.format(RedChisq)))
            PvI_plot.text(0.92, 0.69, Coefficients, verticalalignment='bottom', horizontalalignment='right')
        if weighted=='true':
            Coefficients = '\n'.join((r'$\alpha = {0:.2f} \pm {1:.2f}$'.format(PvIpow[0],RoundErrorPvI),
                                    r'$A = {0:.1f} \pm {1:.1f}$'.format(PvIcoeff,RoundErrorCoeff)))
            PvI_plot.text(0.92, 0.75, Coefficients, verticalalignment='bottom', horizontalalignment='right')

    return PvI_plot

# =============================================================================
# Polarization fraction as a function of estimated total intensity
# =============================================================================
def PvIest(CatPI, CatdPI, CatIest, Iscale=None, Pscale=None, showfit='false', weighted='true', errorbars='false', FigSize=[5.12,3.84], Source=None):  
    # Creating new polarization fractions
    CatP = 100.0 * CatPI / CatIest
    CatdP = CatP*(CatdPI/CatPI)#**2.0
                  
    PvI_plot = plt.figure(figsize=(FigSize[0],FigSize[1])) # Creating the figure object
    plt.xlabel('Stokes Total Intensity $I$ (mJy arcsec$^{-2}$)')
    plt.ylabel("Polarization Fraction $P$ (%)")
    # Plots with and without errorbars 
    if errorbars=='true':
        markers, caps, bars = plt.errorbar(CatIest, CatP, yerr=CatdP, fmt ='.', color='black', ecolor='gray', markersize=1, mfc='none',elinewidth=0.9)
        [bar.set_alpha(0.25) for bar in bars]
    else:
        plt.plot(CatIest, CatP, '.', color='black', markersize=1, mfc='none') 
    plt.xscale('log')
    axes = plt.gca() # Calling the axes object of the figure
    axes.xaxis.set_ticks_position('both') # Adding ticks to each side
    axes.yaxis.set_ticks_position('both') # Adding ticks to each side
    plt.tight_layout() # Using all available space in the plot window
    
    # Checking the automated limits for I and P
    xmin, xmax = plt.xlim()
    ymin, ymax = plt.ylim()
    # Setting mininum and maximum values of the I plot
    if Iscale==None:
        Iscale = np.empty(2)
        Iscale[0] = xmin
        Iscale[1] = xmax
    if Pscale==None:
        Pscale = np.empty(2)
        Pscale[0] = ymin
        Pscale[1] = ymax
    # Setting the minimum and maximum values
    plt.xlim(Iscale[0],Iscale[1])
    plt.ylim(Pscale[0],Pscale[1])


    if showfit=='true':
        # Create the array with ln(I), ln(P), and d(ln(P))
        # The numpy package uses log(x) as the natural logarithm ln
        length = CatP.size
        PvIarray = np.empty([length, 3])
        PvIarray[:,0]=np.log(CatIest)
        PvIarray[:,1]=np.log(CatP)
        PvIarray[:,2]=CatdP/CatP
        # Defining linear function to be fitted by curve_fit
        def lnPvlnI(x, a, b):
            return a*x + b
        # Fitting power law to data
        if weighted=='true':
            PvIpow, PvIcov = optimize.curve_fit(lnPvlnI, PvIarray[:,0], PvIarray[:,1], sigma=PvIarray[:,2])
        else:
            PvIpow, PvIcov = optimize.curve_fit(lnPvlnI, PvIarray[:,0], PvIarray[:,1])
        # Uncertainties on the fit
        PvIerr = np.sqrt(np.diag(PvIcov))
        # Converting the coefficient back to its non-logarithmic value
        PvIcoeff = math.exp(PvIpow[1])
        PvIdcoeff = math.exp(PvIpow[1]) * PvIerr[1]
        # Printing the results of the fit
        print('Power law index: '+str(PvIpow[0])+' ± '+str(PvIerr[0]))
        print('Coefficient: '+str(PvIcoeff)+' ± '+str(PvIdcoeff))

        if Source != None:
            PvI_plot.text(0.92, 0.85, Source, verticalalignment='bottom', horizontalalignment='right', fontweight='bold')

    # Adding power law fit
    if showfit=='true':
        def PvIfunc(x, A, B):
            return A*x**B
            
        # Creating x values
        PvIfuncX = np.arange(Iscale[0],Iscale[1], 0.1)
        # Calculating the y values
        PvIfuncY = PvIfunc(PvIfuncX,PvIcoeff,PvIpow[0])
        plt.plot(PvIfuncX, PvIfuncY, 'red')
        # Rounding up the fit error
        RoundErrorPvI = round_decimals_up(PvIerr[0],2)
        RoundErrorCoeff = round_decimals_up(PvIdcoeff,1)

        # Adding the fit results as a label
        Coefficients = '\n'.join((r'$\alpha = {0:.2f} \pm {1:.2f}$'.format(PvIpow[0],RoundErrorPvI),
                                 r'$A = {0:.1f} \pm {1:.1f}$'.format(PvIcoeff,RoundErrorCoeff)))

        PvI_plot.text(0.92, 0.75, Coefficients, verticalalignment='bottom', horizontalalignment='right')

    return PvI_plot

# =============================================================================
# Polarization fraction as a function of column density
# =============================================================================
def PvN(CatP, CatdP, CatN, Nscale=None, Pscale=None, showfit='false', weighted='true', errorbars='false', FigSize=[5.12,3.84], Source=None):  
    PvN_plot = plt.figure(figsize=(FigSize[0],FigSize[1])) # Creating the figure object
    plt.xlabel('Column Density $N_{H_2}$ (cm$^{-2}$)')
    plt.ylabel('Polarization Fraction $P$ (%)')
    # Plots with and without errorbars 
    if errorbars=='true':
        markers, caps, bars = plt.errorbar(CatN, CatP, yerr=CatdP, fmt ='.', color='black', ecolor='gray', markersize=1, mfc='none',elinewidth=0.9)
        [bar.set_alpha(0.25) for bar in bars]
    else:
        plt.plot(CatN, CatP, '.', color='black', markersize=1, mfc='none')   
    plt.xscale('log')
    axes = plt.gca() # Calling the axes object of the figure
    axes.xaxis.set_ticks_position('both') # Adding ticks to each side
    axes.yaxis.set_ticks_position('both') # Adding ticks to each side
    plt.tight_layout() # Using all available space in the plot window

    # Checking the automated limits for I and P
    xmin, xmax = plt.xlim()
    ymin, ymax = plt.ylim()
    # Setting mininum and maximum values of the I plot
    if Nscale==None:
        Nscale = np.empty(2)
        Nscale[0] = xmin
        Nscale[1] = xmax
    if Pscale==None:
        Pscale = np.empty(2)
        Pscale[0] = ymin
        Pscale[1] = ymax
    # Setting the minimum and maximum values
    plt.xlim(Nscale[0],Nscale[1])
    plt.ylim(Pscale[0],Pscale[1])


    if showfit=='true':
        # Create the array with ln(I), ln(P), and d(ln(P))
        # The numpy package uses log(x) as the natural logarithm ln
        length = CatP.size
        PvNarray = np.empty([length, 3])
        PvNarray[:,0]=np.log(CatN)
        PvNarray[:,1]=np.log(CatP)
        PvNarray[:,2]=CatdP/CatP
        # Defining linear function to be fitted by curve_fit
        def lnPvlnN(x, a, b):
            return a*x + b
        # Fitting power law to data
        if weighted=='true':
            PvNpow, PvNcov = optimize.curve_fit(lnPvlnN, PvNarray[:,0], PvNarray[:,1], sigma=PvNarray[:,2])
        else:
            PvNpow, PvNcov = optimize.curve_fit(lnPvlnN, PvNarray[:,0], PvNarray[:,1])
        # Uncertainties on the fit
        PvNerr = np.sqrt(np.diag(PvNcov))
        # Converting the coefficient back to its non-logarithmic value
        PvNcoeff = math.exp(PvNpow[1])
        PvNdcoeff = math.exp(PvNpow[1]) * PvNerr[1]
        # Printing the results of the fit
        print('Power law index: '+str(PvNpow[0])+' ± '+str(PvNerr[0]))
        print('Coefficient: '+str(PvNcoeff)+' ± '+str(PvNdcoeff))

        if Source != None:
            PvN_plot.text(0.92, 0.85, Source, verticalalignment='bottom', horizontalalignment='right', fontweight='bold')

    # Adding power law fit
    if showfit=='true':
        def PvNfunc(x, A, B):
            return A*x**B
            
        # Creating x values
        PvNfuncX = np.arange(Nscale[0], Nscale[1], 1e20)
        # Calculating the y values
        PvNfuncY = PvNfunc(PvNfuncX,PvNcoeff,PvNpow[0])
        plt.plot(PvNfuncX, PvNfuncY, 'red')
        # Rounding up the fit error
        RoundErrorPvN = round_decimals_up(PvNerr[0],2)
        RoundErrorCoeff = round_decimals_up(PvNdcoeff/1e20,1)
        RoundCoefft = (PvNcoeff/1e20)

        # Adding the fit results as a label
        Coefficients = '\n'.join((r'$\gamma = {0:.2f} \pm {1:.2f}$'.format(PvNpow[0],RoundErrorPvN),
                                 ' '))

        PvN_plot.text(0.92, 0.75, Coefficients, verticalalignment='bottom', horizontalalignment='right')

    return PvN_plot

# =============================================================================
# Polarization fraction as a function of temperature
# =============================================================================
def PvT(CatP, CatT, CatdP=None, Tscale=None, Pscale=None, errorbars='false', FigSize=[5.12,3.84], Source=None, Position='right'):  
    PvT_plot = plt.figure(figsize=(FigSize[0],FigSize[1])) # Creating the figure object
    plt.xlabel('Dust Temperature $T_d$ (K)')
    plt.ylabel('Polarization Fraction $P$ (%)')
    # Plots with and without errorbars 
    if errorbars=='true':
        markers, caps, bars = plt.errorbar(CatT, CatP, yerr=CatdP, fmt ='.', color='black', ecolor='gray', markersize=1, mfc='none',elinewidth=0.9)
        [bar.set_alpha(0.25) for bar in bars]
    else:
        plt.plot(CatT, CatP, '.', color='black', markersize=1, mfc='none')   
    plt.plot(CatT, CatP, '.', color='black', markersize=1, mfc='none')
    #plt.xscale('log')
    axes = plt.gca() # Calling the axes object of the figure
    axes.xaxis.set_ticks_position('both') # Adding ticks to each side
    axes.yaxis.set_ticks_position('both') # Adding ticks to each side
    plt.tight_layout() # Using all available space in the plot window

    # Checking the automated limits for I and P
    xmin, xmax = plt.xlim()
    ymin, ymax = plt.ylim()
    # Setting mininum and maximum values of the I plot
    if Tscale==None:
        Tscale = np.empty(2)
        Tscale[0] = xmin
        Tscale[1] = xmax
    if Pscale==None:
        Pscale = np.empty(2)
        Pscale[0] = ymin
        Pscale[1] = ymax
    # Setting the minimum and maximum values
    plt.xlim(Tscale[0],Tscale[1])
    plt.ylim(Pscale[0],Pscale[1])

    if Source != None:
        if Position == 'right':
            PvT_plot.text(0.92, 0.85, Source, verticalalignment='bottom', horizontalalignment='right', fontweight='bold')
        if Position == 'left':
            PvT_plot.text(0.18, 0.85, Source, verticalalignment='bottom', horizontalalignment='left', fontweight='bold')

    return PvT_plot

# =============================================================================
# Comparison plot between measured and predicted fluxes at 214 µm
# =============================================================================
def IvIest(CatI, CatdI, CatIest, Iscale=[0,65], FigSize=[5.12,3.84], Source=None, Position='top'):

    # Creating reference line
    linemin = Iscale[0]
    linemax = Iscale[1]
    RefX = np.arange(linemin, linemax, 0.1)
    RefY = np.arange(linemin, linemax, 0.1)

    # Fitting a linear relation to IvI
    def IvIfunc(x, a, b):
        return a*x**b
    IvIpara, IvIcov = optimize.curve_fit(IvIfunc, CatIest, CatI)#, sigma=CatdI)
    Coeff_1000 = IvIpara[0]*1000

    print('Coefficient: '+str(Coeff_1000)+r' $\times 10^{-3}$')
    print('Exponent: '+str(IvIpara[1]))
            
    # Creating x values
    IvIfuncX = np.arange(linemin, linemax, 0.1)
    # Calculating the y values
    IvIfuncY = IvIfunc(IvIfuncX,IvIpara[0],IvIpara[1])

    # Fitting a linear relation to IvI
    def IvIfunc2(x, b):
        return x + b
    IvIpara2, IvIcov2 = optimize.curve_fit(IvIfunc2, CatIest, CatI)#, sigma=CatdI)

    print('Constant: '+str(IvIpara2[0]))
            
    # Creating x values
    IvIfuncX2 = np.arange(linemin, linemax, 0.1)
    # Calculating the y values
    IvIfuncY2 = IvIfunc2(IvIfuncX2,IvIpara2[0])


    # IvI
    PvI_plot2 = plt.figure(figsize=(FigSize[0],FigSize[1])) # Creating the figure object
    plt.xlabel('Predicted 214 $\mu$m Flux Density (mJy arcsec$^{-2}$)')
    plt.ylabel('Measured 214 $\mu$m Flux Density (mJy arcsec$^{-2}$)')
    plt.plot(CatIest, CatI, '.', color='black', markersize=1, mfc='none')
    plt.plot(RefX, RefY, 'orange')
    plt.plot(IvIfuncX, IvIfuncY, 'red')
    plt.plot(IvIfuncX2, IvIfuncY2, color='green', linestyle='dashed')
    #plt.xscale('log')
    #plt.yscale('log')
    axes2 = plt.gca() # Calling the axes object of the figure
    axes2.xaxis.set_ticks_position('both') # Adding ticks to each side
    axes2.yaxis.set_ticks_position('both') # Adding ticks to each side
    axes2.set_xlim(linemin, linemax)
    axes2.set_ylim(linemin, linemax)
    plt.tight_layout() # Using all available space in the plot window

    # Positioning legend
    if Position == 'top':
        xpos = 0.2
        ypos1 = 0.85
        ypos2 = 0.7
        align = 'left'
    if Position == 'bottom':
        xpos = 0.92
        ypos1 = 0.35
        ypos2 = 0.2
        align = 'right'

    if Source != None:
            PvI_plot2.text(xpos, ypos1, Source, verticalalignment='bottom', horizontalalignment=align, fontweight='bold')
    # Adding the fit results as a label
    Coefficients = '\n'.join((r'$A_1 = {0:.2f}$'.format(Coeff_1000)+r' $\times 10^{-3}$',
                              r'$\alpha_1 = {0:.2f}$'.format(IvIpara[1]),
                              r'$C_2 = {0:.2f}$'.format(IvIpara2[0]),))

    PvI_plot2.text(xpos, ypos2, Coefficients, verticalalignment='bottom', horizontalalignment=align)

    return PvI_plot2

# =============================================================================
# Function to create an Histogram and fit a Gaussian function
# =============================================================================
def Histogram(CatdO, binsize=5.0, showfit='no', FigSize=[5.12,3.0], Source=None, LabelSide='left', FilO=None, PlanckO=None):
    # Size of the catalog
    CatSize = CatdO.shape[0]
    
    # Method to create an histogram from a vector catalog
    print('Calculating mean and standard deviation of the catalog')
    print()
    
    # Conversion from degrees to radians
    conv_rad = math.pi/180.0
    conv_deg = conv_rad**-1.0
    
    # Calculating the regular mean
    norm_mean = np.mean(CatdO)
    norm_std = np.std(CatdO)
    
    print('Mean: ' + str(norm_mean))
    print('Standard deviation: ' + str(norm_std))
    print()
    
    # Calculating the circular mean and circular standard deviation
    # assuming boundaries of -90 to 90 degrees
    circ_mean = conv_deg * stats.circmean(conv_rad*CatdO, 
                                            low=-math.pi/2.0, 
                                            high=math.pi/2.0)
    circ_std = conv_deg * stats.circstd(conv_rad*CatdO, 
                                        low=-math.pi/2.0,  
                                        high=math.pi/2.0)
    
    print('Circular mean: ' + str(circ_mean))
    print('Circular standard deviation: ' + str(circ_std))
    print()
    
    # Calculating the number of bins by rounding to lowest integer
    number_bins = math.floor(180.0 / (1.0 * binsize))
    eff_binsize = 180.0/number_bins
    
    print('Number of bins: ' + str(number_bins))
    print('Bin size used: ' + str(eff_binsize))
    print()
    
    # Creating the histogram arrays
    histo_bins = np.zeros(number_bins)
    histo_values = np.zeros(number_bins)
    # Identifying the values of the bins
    for i in range (0,number_bins):
        histo_bins[i] = (i + 0.5) * eff_binsize - 90.0
    # Counting the number of vectors per bin - BRUTE FORCE VERSION
    # CONSIDER USING NUMPY HISTOGRAM FUNCTION INSTEAD
    for i in range (0,number_bins):
        # Limits of the bin
        min_limit = histo_bins[i] - 0.5 * eff_binsize
        max_limit = histo_bins[i] + 0.5 * eff_binsize
        # Checking every vector if they should be in the bin
        for j in range (0,CatSize):
            if (CatdO[j] >= min_limit) and (CatdO[j] < max_limit):
                histo_values[i] = histo_values[i] + 1
    
    # Find array value closest to circular mean    
    abs_array = np.abs(histo_bins - circ_mean)
    smallest_diff = abs_array.argmin()        
    # Calculating shift needed to center the histogram on circular mean
    shift = math.ceil(number_bins/2.0) - smallest_diff
    # Rolling the histogram
    print('Using the circular mean to center the histogram')
    print('Rolling histogram by '+str(shift)+' elements')
    print()
    histo_bins_rolled = np.roll(histo_bins, shift)
    histo_values_rolled = np.roll(histo_values, shift)
    
    # Fixing values of rolled bins so array stays sorted
    if shift > 0:
        for i in range(0, shift):
            histo_bins_rolled[i] = histo_bins_rolled[i] - 180.0
    if shift < 0:
        for i in range(number_bins+shift, number_bins):
            histo_bins_rolled[i] = histo_bins_rolled[i] + 180.0
    
    # Attempting Gaussian fit to the data
    # Defining the Gaussian function, stolen from the internet
    # https://lmfit.github.io/lmfit-py/model.html
    def gaussian(x, amp, cen, wid):
        y = amp*np.exp(-(x-cen)**2 / (2*wid**2))
        return y
    
    print('Fitting a Gaussian profile to the histogram')
    print()
    
    # Calling the lmfit package to model the 'gaussian' function                        
    gaussian_model = Model(gaussian)
    gaussian_fit = gaussian_model.fit(histo_values_rolled, 
                                        x=histo_bins_rolled, 
                                        amp=np.mean(histo_values_rolled), 
                                        cen=circ_mean, wid=circ_std)
    # Creating the array containing the parameters
    fit_params = np.empty([3,2])
    # Amplitude
    fit_params[0,0] = gaussian_fit.params['amp'].value
    fit_params[0,1] = gaussian_fit.params['amp'].stderr
    # Mean
    fit_params[1,0] = gaussian_fit.params['cen'].value
    fit_params[1,1] = gaussian_fit.params['cen'].stderr
    # Standard deviation
    fit_params[2,0] = gaussian_fit.params['wid'].value
    fit_params[2,1] = gaussian_fit.params['wid'].stderr

    print('Goodness of fit')
    print('Number of iterations: '+str(gaussian_fit.nfev))
    print('Reduced chi-squared: '+str(gaussian_fit.redchi)) 
    print()
    
    print('Fitted parameters')
    print('Amplitude:          '+str(gaussian_fit.params['amp'].value))
    print('                  ± '+str(gaussian_fit.params['amp'].stderr))
    print()
    print('Mean:               '+str(gaussian_fit.params['cen'].value))
    print('                  ± '+str(gaussian_fit.params['cen'].stderr))
    print()
    print('Standard deviation: '+str(gaussian_fit.params['wid'].value))
    print('                  ± '+str(gaussian_fit.params['wid'].stderr))
    print()
    
    # Plotting the histogram
    histo_plot = plt.figure(figsize=(FigSize[0],FigSize[1])) # Creating the figure object
    plt.bar(histo_bins_rolled, histo_values_rolled, 
                    width=eff_binsize, fill=False) # Creating the bar plot
    plt.xlabel(r'Polarization Angle $\theta$ (Degree)')
    plt.ylabel('Number')
    plt.tight_layout() # Using all available space in the plot window
    if Source != None:
        # Creating the reference box
        props = dict(boxstyle='square', facecolor='white', alpha=0.8) # Setting mild transparency
        # Including the information
        textstr = '\n'.join((
            Source,
            r'$\overline{\theta}=%.1f^{\circ}$' % circ_mean,
            r'$\sigma_{\theta}=%.1f^{\circ}$'% circ_std))
        if FilO != None:
            textstr = '\n'.join((
            Source,
            r'$\overline{\theta}=%.1f^{\circ}$' % circ_mean,
            r'$\sigma_{\theta}=%.1f^{\circ}$'% circ_std,
            r'$\phi_{F}=%.1f^{\circ}$'% FilO))
        if PlanckO != None:
            textstr = '\n'.join((
            Source,
            r'$\overline{\theta}=%.1f^{\circ}$' % circ_mean,
            r'$\sigma_{\theta}=%.1f^{\circ}$'% circ_std,
            r'$\overline{\theta}_{P}=%.1f^{\circ}$'% PlanckO))
        if FilO != None and PlanckO !=None:
            textstr = '\n'.join((
            Source,
            r'$\overline{\theta}=%.1f^{\circ}$' % circ_mean,
            r'$\sigma_{\theta}=%.1f^{\circ}$'% circ_std,
            r'$\overline{\theta}_{P}=%.1f^{\circ}$'% PlanckO,
            r'$\phi_{F}=%.1f^{\circ}$'% FilO))
        if LabelSide == 'left':
            histo_plot.text(0.17, 0.75, textstr, verticalalignment='center', horizontalalignment='left', bbox=props, fontweight='bold')
        if LabelSide == 'right':
            histo_plot.text(0.94, 0.75, textstr, verticalalignment='center', horizontalalignment='right', bbox=props, fontweight='bold')
    # Setting range
    xmin = histo_bins_rolled[0]-eff_binsize/2.0
    xmax = histo_bins_rolled[number_bins-1]+eff_binsize/2.0
    plt.xlim(xmin, xmax)
    # Modifying the axes directly
    axes = plt.gca() # Calling the axes object of the figure
    axes.xaxis.set_ticks_position('both') # Adding ticks to each side
    axes.yaxis.set_ticks_position('both') # Adding ticks to each side
    # Plotting the best fit Gaussian
    # Creating x values
    gaussian_bins = np.arange(xmin, xmax, 0.1)
    # Calculating the y values
    gaussian_values = gaussian(gaussian_bins, fit_params[0,0], 
                                fit_params[1,0], fit_params[2,0])
    # Plotting the Gaussian fit
    if showfit == 'yes':
        plt.plot(gaussian_bins, gaussian_values, 'k')

    # Adding reference lines
    plt.axvline(x=-180.0, color='black', linestyle='--', linewidth=2)
    plt.axvline(x=-90.0, color='black', linestyle='--', linewidth=2)
    plt.axvline(x=0.0, color='black', linestyle='--', linewidth=2)
    plt.axvline(x=90.0, color='black', linestyle='--', linewidth=2)
    plt.axvline(x=180.0, color='black', linestyle='--', linewidth=2)
    plt.axvline(x=circ_mean, color='red', linestyle='-', linewidth=2)
    if FilO != None: 
        plt.axvline(x=FilO, color='blue', linestyle='-.', linewidth=2)
    if PlanckO != None: 
        plt.axvline(x=PlanckO, color='green', linestyle='-.', linewidth=2)
    
    # Returning figure and fits parameters
    return histo_plot, fit_params

# Function to round error decimals up
#  See https://kodify.net/python/math/round-decimals/
def round_decimals_up(number:float, decimals:int=2):
    """
    Returns a value rounded up to a specific number of decimal places.
    """
    if not isinstance(decimals, int):
        raise TypeError("decimal places must be an integer")
    elif decimals < 0:
        raise ValueError("decimal places has to be 0 or more")
    elif decimals == 0:
        return math.ceil(number)

    factor = 10 ** decimals
    return math.ceil(number * factor) / factor